data {
  
  int n_items;       // number of items
  int n_respondents; // number of respondents
  int n_obs;         // number of observations
  
  int k[n_obs];      // respondent ID
  int j[n_obs];      // item ID
  real y[n_obs];     // item placement by respondent
  
}

parameters {
  
  vector[n_respondents] alpha_raw;
  real alpha_mu;
  real<lower = 0> alpha_sigma;
  
  vector[n_respondents] beta_raw;
  real<lower = 0> beta_mu;
  real<lower = 0> beta_sigma;
  
  vector[n_items-2] zeta_raw;
  
  real<lower = 0> sigma[n_items];
  
}

transformed parameters {
  
  vector[n_items] zeta;
  vector[n_respondents] alpha;
  vector[n_respondents] beta;
  
  zeta[1] = -1;
  zeta[2] = 1;
  zeta[3:n_items] = zeta_raw;

  // Non-centered parameterization for efficiency
  for(m in 1:n_respondents) {
    alpha[m] = alpha_mu + alpha_raw[m] * alpha_sigma;
    beta[m] = beta_mu + beta_raw[m] * beta_sigma;
  }
  
}

model {
  
  vector[n_obs] y_star;
  vector[n_obs] sigma_star;
  // int pos;
  
  // Priors
  alpha_mu ~ normal(0, 5);
  alpha_sigma ~ cauchy(0, 2.5);
  alpha_raw ~ normal(0, 1);
  
  beta_mu ~ normal(0, 5);
  beta_sigma ~ cauchy(0, 2.5);
  beta_raw ~ normal(0, 1);
  
  zeta_raw ~ normal(0, 5);

  sigma ~ cauchy(0, 2.5);
  
  // Likelihood
  
  // pos = 1;
  for(m in 1:n_obs) {
    if(j[m] != n_items+1) {
      y_star[m] = alpha[k[m]] + beta[k[m]] * zeta[j[m]];
      sigma_star[m] = sigma[j[m]];
    }
  }
  
  y ~ normal(y_star, sigma_star);

}

generated quantities {
}
